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Bayesian approach to gravitational lens model selection: 
constraining H with a selected sample of strong lenses 
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ABSTRACT 

Bayesian model selection methods provide a self-consistent probabilistic framework to test 
the validity of competing scenarios given a set of data. We present a case study application 
to strong gravitational lens parametric models. Our goal is to select a homogeneous lens sub- 
sample suitable for cosmological parameter inference. To this end we apply a Bayes factor 
analysis to a synthetic catalog of 500 lenses with power-law potential and external shear. For 
simplicity we focus on double-image lenses (the largest fraction of lens in the simulated sam- 
ple) and select a subsample for which astrometry and time-delays provide strong evidence for 
a simple power-law model description. Through a likelihood analysis we recover the input 
value of the Hubble constant to within 3er statistical uncertainty. We apply this methodology 
to a sample of double image lensed quasars. In the case of B1600+434, SBS 1520+530 and 
SDSS J1650+4251 the Bayes' factor analysis favors a simple power-law model description 
with high statistical significance. Assuming a flat ACDM cosmology, the combined likeli- 
hood data analysis of such systems gives the Hubble constant Hq = 76I5 5 kms~ 1 Mpc _1 
having marginalized over the lens model parameters, the cosmic matter density and consis- 
tently propagated the observational errors on the angular position of the images. The next 
generation of cosmic structure surveys will provide larger lens datasets and the method de- 
scribed here can be particularly useful to select homogeneous lens subsamples adapted to 
perform unbiased cosmological parameter inference. 
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1 INTRODUCTION 

Strong gravitational lenses are powerful cosmological probes. 
The distorted images of lensed background sources carry infor- 
mation on the distribution of invisible matter in cosmic struc- 
tures. Furthermore, as light-rays from distant variable sources 
propagate along differently distorted paths, time-delays be- 
tween luminosity variations of the source images test the un- 
derlying cosmological expansio n (for exhaustive review see 
Schnei der. Ehlers & Falco 1991 ; Petters, Levine & Wamb sgansj 
200ll : ls"chneider. Kochanek & Wam bsganss 2006). 

The use of le ns time-delays as cos mic standard rulers was ini- 
tially proposed bv lRefsdall d 1964 1 19661) . To date time-delays have 
been measured only in 21 lensed quasars out of the few hundreds 
known strong lens systems. Despite the scarcity of observations, 
lots of effort has been devoted to using these data to infer the value 
of the Hubble constant, Ho - This is because time-delays, differently 
from other standard cosmological tests, do not rely on local cali- 
bration measurements, such as the distanc e ladder method (for an 
alternative method see lChavez et alj|2012l) . In the next decade, ad- 
ditional lens time-delay data will become available thanks to obser- 
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vational programs such as the COSMOGRAIlJj project collabora- 
tion, as well as the ILMT projecQ Furthermore, the next generation 
of cosmic structure surveys such as the Dark Energy Survey (DES), 
the Large Synoptic Telescope (LSST) or the EUCLID satellite mis- 
sion will detect large sample of strong lens systems for which 
time-delays can be accurately measured through follow-up obser- 
vations. Time-delays from such datasets will be particularly useful 
to infer cosmological constraints that are compl ementary to those 
obtained from other cosmic probes (see e.g. iDobke et al.l 120091 : 
ICoe & Moustakasll2009l : IOguri & MarshalfeOloHLindeifeOllI) . 

Time-delay measurements are particularly challenging be- 
cause they require an intense monitoring of the lens system. How- 
ever, the main limitation to fully exploit the cosmological infor- 
mation encoded in such data arises from the uncertainty on the lens 
mass distribution as well as the presence of perturbing masses along 
the line-of-sight. Because of t his, several analyses have foc used on 
"golden lenses" (Press 1996; Wi lliams & Schechterlll997l) . These 
are systems for which there is a sufficient number of observational 
features (e.g. presence of arcs, multiple point-like images, flux ra- 
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tios measurements) such as to constrai n the lens potential indepen- 
dently of time-delays. For example, Wucknitz, Bigg s & Browne! 
(2004) have studied the lens system B0218+357 and inferred 
H = 78 ± ekms^Mpc" 1 for a fiat Cold Dark Matter model 
with Cosm ological Constant (ACDM) with mean matter density 
n m = 0.3. ISuvuetai]j2010h have analysed the lens B 1608+656 
and for the same cosmology they have obtained Ho = 70.6 ± 
3.1kms _1 Mpc _1 . 

Alternatively one can infer cosmological parameter con- 
straints u sing a statistical s a mple of lens time-delay measurements 
(see e.g. ISaha et al.l 120061 : lOguril 120071) . Such an approach still 
requires modeling the gravitational potential of each lens in the 
sample. Nevertheless, one can hope that systematics due to in- 
dividual mass model uncertainties are averaged out. In such a 
case sample selection effects can be the main source of error 
dOguri. Keeton & Dalai 2051 ; I Ogurill 20061) . 

Modeling the lens mass d i stribut i on can be distinguished 
in parametric (e. g. see lOguril 120021 : iKeeton. Gaudi & Pettersl 



l2005h . Bayesian model selection methods have already been ap- 
pli ed to a variety of problems in cosmology and astrophysics (e.) 
see|jaffdl996l:lMarshall. Hobson & Slosarl2003l:ISaini et al.l2004 
Bassett, Corasaniti & Kunz 2008; Mukherjee, Parkinson & Liddle 



2006; Mukherjee et al. 200' 



Trottal 



2007 



Ford & Gregory 



20031: lQgurill2007l) and non-paramet r ic (e.g . see Kochanek 1991; 
Saha & Williams 1997J; iKoopmansl 120051 ; IVegetti & Koopmana 
20091 : ISuvu etal]l2009h methods. The latter uses linear inversion 



algorithms to constrain the lens potential directly from the inten- 
sity images of the lens system, while the former uses the measured 
properties of the lens to const rain a parametrized form of the lens 
potential. A third approach by Alard (2007, 2008) reconstructs the 
lens potential as a perturbative expansion around the Einstein ring 
of the lens system. 

Non-parametric methods have been extensively used in a vast 
literature. As an example a non-p arametric reconstruction algo- 
rithm is implemented in Pixeleno Williams & S aha 2000), a nu- 
merical code com monly used for lens studies. Using this code 
ISaha et al.l d2006l) have found H = 72±l 1 kms^Mpc -1 from 
a sample of 10 time-delay lenses f or a flat ACDM with Q m = 0.3. 
Recently Paraficz & Hjorfrj d2010h have performed the same anal- 
ysis on an extended sample of 18 lenses and obtained Ho = 

66^4 kms _1 Mpc _1 . A statistical analysis of a sample of lensed 2 GRAVITATIONAL LENSES 



Cornish & Littenberd 120071; IGregorv & Fischerll2010l) . In the lit 
erature, the analysis of Bayes factors has been applied to 
non-parametric lens recons truction techniques dSuvuetal.ll2006l : 
IVegetti & Koopmansl [2009 ) as well a s selecting comp licate para- 
metric lens models of galaxy clusters djullo et al.ll2007l) . 

Here, we present a case study of the use of Bayes factors to 
construct a homogeneous sample of double lensed quasars to de- 
rive bounds on Ho - By homogeneous we mean a sample consisting 
of lenses whose data, astrometry and time-delays of the images, 
provide evidence for the same lens model description at the same 
level of statistical significance. Differently from previous statisti- 
cal approaches we consistently propagate all observational uncer- 
tainties on the posterior probabilities, including errors on the angu- 
lar position of the images, as well as marginalizing over nuissance 
model parameters rather than assuming hard priors. The possibil- 
ity of selecting homogeneous lens data to infer cosmological pa- 
rameter constraints will become particularly important for future 
survey programs which will detect a large number of strong grav- 
itational lenses. We argue that the use of Bayes factors provides 
a self-consistent probabilistic method to build subsample of data 
which are not dominated by astrophysical selection effects. 

The paper is organized as follows: in Section [2] we review the 
lens equations and the modeling of the gravitational lens potential, 
while in Section [3] we discuss the Bayesian model selection. We 
present the results of the lens model selection analysis on simulated 
data in Section [4] and on real data in Section [5] In Section [6] we 
describe the results of the cosmological parameter inference and 
present our conclusions in Section|7] 



quasars using a parametric approach has been performed by a num- 
ber of authors. For instanc e, assuming an isothermal lens potential 
iGiovi & Amendok] ( EoOlh have inferred li mits on Hp for differ- 
ent cosmological models. The analysis bv lOguril d2007l) has prop- 
agated the uncertainties in the lens model parameters and found 
Ho = 70 ± 6 km s _1 Mpc _1 from a sample of 16 lenses for a flat 
ACDM with Q. m = 0.26. Whether using non-parametric lens re- 
construction algorithms or a parametric model a key aspect of these 
analyses is the assessment of what constitutes a good description of 
the lens potential given the data. 

Independently of the approach used, distinguishing between 
the potentially infinite number of possibilities has been mostly 
based on x 2 -statistics. However, parameter fitting only establishes 
how well a model reproduce the data for a given set of model pa- 
rameter values. Deciding whether one model is preferable over an- 
other is a question of model selection rather than quality of pa- 
rameter fit. In other words, which model has a higher probability 
of being the correct model description of the observations? Do the 
data justify a more complex description of the system (additional 
parameters)? 

In the Bayesian statistical framework these questions are ad- 
dressed by the analysis of the Bayesian evidence and the eval- 
uation of Bayes factor ( [Jeffreys! Il96ll ; iMacKavl 120031 ; IGregorv! 



www.qgd.uzh.ch/projects/pixelens/ 



2.1 Lens equations 

Here, we briefly review the basic equations describing the forma- 
tion of images in strong gravitational lenses. We assume that the 
source is lensed by an object that can be treated as a single lens 
plane. 

Let us consider an angular coordinate system centered on a 
lens at redshift z\ and a source with angular position (3 at red- 
shift z s . According to Fermat's principle light-rays from the source 
follow paths that extremize the arrival-time. Let t(G, f3) be the 
arrival-time of rays observed at an angle 0. We can estimate this 
function by considering the geodesies of photons connecting the 
source-lens and the lens-observer planes respectively. Two effects 
contribute to the arrival-time, a geometrical term which accounts 
for the different length of the lensed paths and a gravitational 
term due to the Shapiro ef fect (for a detailed derivation e.g. see 
iBlandford & Naravanll 986). The geometrical contribution reads as 

1 + zi Din 



a) 



2 cAs 

where c is the speed of light, while Di, D s and Di s are the angular 
diameter distances between observer and lens, observer and source, 
and lens and source respectively; in the following, we note 6 and /3 
the norms of and f3. The gravitational term is given by 



{0,(3) 



(2) 
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where is a 2-D gravitational potential generated by the pro- 

jected surface mass density on the lens plane at the angular position 
of the source, 2(0), as given by the Poisson-like equation 



V 2 ¥(0) = GE(0), 



(3) 



where G is the Newton gravitational constant. It is useful to in- 
troduce the critical density E c = p p . tne convergence 
«: = E/S c and the reduced potential ij) = 2¥/(GS c ) to write 
the total arrival-time as 



t(0,p) 



(1 + Z1 ) 



AA 

cDis 



/3) 2 -V(0) 



(4) 



from which we can derive all relevant equations describing the 
properties of the images. 

Lens Equation: let us assume the source to be fixed at the 
angular location p, the arrival-time only depends on and images 
will form at extrema of the arrival time corresponding to \jt — 0. 
Using Eq. {4} we obtain the lens equation 



(3 = 6- V ^(0) 



(5) 



Magnification: gravitational lensing deforms the images of 
a lensed source and since the luminosity per surface area is con- 
served, images can appear magnified or demagnified. The ampli- 
tude of the effect can be quantified by considering the relation be- 
tween the coordinates on the source plane and the lens one. Let 



A = yp=X-yyij}, 



(6) 



be the Jacobian matrix of the angular coordinate transformation, 
this relation implies that a solid-angle element S/3 2 of the source is 
mapped to a solid-angle element of the image S6 2 by the inverse of 
the determinant of XjP- Hence, the magnification is given by 



/ ' '■■ 



S/3 2 



det.4' 



(7) 



Time-delay: let us consider a double-image lensed source 
characterized by a time-variable luminosity, the time-delay in the 
appearance of the luminosity variation between two images A and 
B is given by: 

A A 



AtAB = t(O A ,P)-t(0B,P) = (l+Z 1 ) 



cDu 



\{9a 



Pf 



•iK0A)-i(fe 



Pf+^(0B) 



(8) 



The position of the images, their relative flux ratios and time- 
delays are the main observable features of a strong lens systems. 



2.2 Lens models 

Our goal is to select a homogeneous sample of simple lens systems 
that can be easily observed in large sky surveys. For this reason we 
focus on individual galaxies which lens distant quasar sources. 



2.2.7 Power-law density profile 

The mass distribution of lens spiral and elliptical gal axies is wel l 
approximated by power-law density profiles (e.g. see lRusinll2003h 
for which the lens potential assume the form: 



(9) 



where b is a deflection scale. The singular isothermal sphere (SIS) 
model teinnev & Trem aine 1987) corresponds to n = 2, for which 
6 = AivDisO- 2 / D B , where a is the velocity dispersion of the galaxy. 
Measurements of galaxy density profiles indicates that the slope 
parameter n is generally close to the isothermal value, though 
some systems have revealed shallow profiles with n < 1 (e.g. see 
ISalucci et a l. 2007). By computing the second derivative of Eq. {9]l 
with respect to 6 we obtain the convergence: 



<0) = 



2-n 



(10) 



In the case of a double image lens with point-like images located at 
6a and 9b, the deflection scale can be written as: 



b = 



+ I 



+ i 



(11) 



where we have used Eq. (0. For such a system the time-delay be- 
tween the two images is given by l&chanek 20(1): 



AtAB = (1 + Zl) 



cAs 



x (!-(«}) 



+ 



B log ( 



K0)-(«>]iog — )dde 



where 



+ 



(12) 



(13) 



is the mean surface mass density in the ring between 9a and 8b- 
The above formulae can be rewritten using Eq. (QJJ): 



AtAI 



(l + zi) 



AA 

cA„ 



Zi3 — n/i2 



1 / \ 2 " 

2- <K) 3" 



(1 - (k))0 a 0e 



(14) 



and 



9l + e\ 



(15) 



which explicitely depend on the lens model parameters. 

Similarly we can write the flux-ratios between images A and B 
in terms of their angular positions and the parameters of the power 
law potential using Eq. 10. The flux-ratio is given by 



TAB 



Fa 
F B ' 



where 
1 

F 



l-(2 



(16) 



(17) 



with i = A, B. This accounts only for luminosity differences be- 
tween images due to lensing magnification. ^From Eq. dl7| l we can 
see that flux-ratio measurements can strongly constrain the lens 
model parameters, thus in combination with time-delays such data 
can reduce or break the mass-sheet degeneracy and lead to tighter 
bounds on cosmological parameters. 

Several source of systematic errors currently limit the use of 
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flux-ratio measurements. In fact, these can be affected by inhomo- 
geneous dust extinction across the angular size of the lens sys- 
tem, thus altering the contribution of the lens magnification en- 
coded in the flux-ratio measurements. In principle dust effects can 
be taken into account through color analysis. On the other hand, 
rapidly varying lensing effects can be a more important source 
of systematic uncertainty that remains hard to model. As shown 
by ISchild & Smithl il99lh in the case of the double image lens 
Q0957+561 flux-ratios can vary over time. Fine variations on short 
time scales can be caused by microlensing of the structure of the 
lensed quasar, while trends over longer periods are more indicative 
of the mass spectrum of the lens galaxy. 

Radio measurements of flux ratios are primarely affected by 
milli-lensing events, while in the optical both milli-lensing and 
micro-lensing causes flux anomalies. In the latter case flux-ratio 
estimates from spectroscopic measurements may alleviate the sys- 
tematic effect du e to micro-lensing as recently pointed out by 
ISluse etaljj2012h . 

Henceforth, it is not surprising that flux-ratios can be 
an effective probe of the lens m ass distribution (e.g. see 
iGoicoechea. Gil-Merino &"u ilan 2005), but at the same time if un- 
modelled, flux anomalies introduce dominant systematic errors in 
the cosmological parameter inference. Overcoming these limita- 
tions requires a systematic monitoring of the lens systems phased 
on the lens time-delay as well as an accurate modeling of the lens 
inner structure. Such detailed analyses have yet to be performed and 
are beyond the scope of our work. Because of this cosmological 
constraints inferred in combination with flux-ratio measurements 
should be taken with a grain of salt. In Section 1531 we will show 
results obtained by combining time-delays with flux ratios only for 
illustrative purposes. 



Eq. d 1 8b adds two parameters to the lens model. Consequently, 
Eq. dl4t-dl7b are modified by terms which explicitely depends on 
the external shear parameters. Such terms can be computed analyt- 
ically using Eq. dl8b . their derivation is quite cumbersome and we 
leave it to Appendix [A] An important aspect of the presence of ex- 
ternal shear concerns the geometrical configuration of the images. 
In particular the angular separation of the images Gab = 8b — 8 a, 
the asymmetry -Rab = \(8b — 6 a) /(8b + 8a)\, and the the de- 
viation from colinearity e = j @ab — 180|/180 are indicators of 
pert urbations about a monopole-like lens poten tial. As an exam- 
ple, ISchneider. Kochanek & Wambsgansll d2006T) have shown that 
the time-delays have little sensitivity to the quadrupole term if the 
images are opp osite to each other with respect to the lens posi- 
tion. Similarly, |Oguri (2007) has shown that for image pairs with 
an asymmetry parameter Rab > 0.2 the time-delay is not signif- 
icantly influenced by an external shear, a third-order external per- 
turbation or the presence of dark matter subhalos. In contrast, for 
such systems the time-delay is more sensitive to the radial slope of 
the potential (which quantifies the deviation from isothermality), 
which appears to be the case especially for image pairs with a wide 
angular separation. This offer a first empirical criterion to select a 
homogeneous lens sample. 

In addition to external shear, a non-spherical lens mass dis- 
tribution can lead to a non-trivial angular dependence of the lens 
potential that contributes to the time-delay. If unmodelled this also 
introduce a systematic error altering the cosmological parameter 
inference. In the case of lenses with images lying opposite to each 
other the contribution from this internal shear appears onl y as a 
second order correction to the time-delay Kochanek (20021. This 
suggests that in general the presence of internal shear should not be 
neglected. 



2.2.2 External shear 

The presence of mass perturbators outside the lens system can di- 
rectly alter the lens potential and introduce uncertainties in the 
modeling of the lens potential. This is especially the case if a sin- 
gle lens galaxy is not isolated, rather is part of a group or a cluster 
of galaxies where an external shear field can deform the monopole 
potential Eq. {5). Such deformation can be modeled by adding a 
quadrupole term and in polar angular coordinates — (8, <j>) this 
reads as 



-i 7 2 cos2( 



(18) 



where 7 is an amplitude parameter which measure the strength of 
the external shear and 7 its direction. Observational constraints 
on the external shear have been mainly inferred from the study of 
multiple image lenses. As an example, a number of studies of in- 
dividual quadruple lens systems hav e shown that the shear ampli- 
tude is rather large, 7 ~ 0.1 — 0.3 dFischer. Schade & Barrientosl 
1 19981: iKneib. Cohen & Hiortrtood) . In contrast, earlier expecta- 
tions were in the range ~ 0.02 — 0.05 (Keeton, Kochanek & Seljak 
1 1997b . Using results from N-body simulati ons in combination 
with s emi-analytic models of galaxy formation lHolder & Schechteil 
42003b have shown that distribution of values of 7 is nearly- 
Gaussian with a p eak at 7 ~ 0.1 and a rapid decay for 7 > 0.2. 
On the other hand lWong et al.l d201 lh have performed an accurate 
study of nine lens systems in galaxy groups and clusters which in- 
dicates that on average 7 = 0.08, with a distribution ranging from 
0.02 to 0.17. 

Accounting for the effect of external shear as described by 



3 BAYESIAN STATISTICAL ANALYSIS 

Our goal is to use Bayesian model selection to construct a statisti- 
cally homogeneous sample of gravitational lenses from which to in- 
fer cosmological constraints using time-delays. More specifically, 
we aim to select lenses for which data provide strong evidence in 
favor of a simple power-law model description and exclude with 
high statistical significance the presence of external shear which 
might contribute to sample selection effects in single lens galaxies. 
Here we review the basic elements of Bayesian statistical analysis. 



3.1 Model parameter estimation and error propagation 

Let us consider a set of observations of a double image lens D, 
consisting of the angular position of the images 8i (i = A, B), the 
time-delay AtAB (and eventually flux ratio tab). We want to con- 
strain a lens model A4 specified by a set of model parameters ex. 
For the time being let us assume hard priors on the cosmological pa- 
rameters which determine the cosmic angular distances in Eq. (8). 
The first step of the data model comparison consists of inferring 
the best fit values of the lens model parameters as well as their un- 
certainties expressed in terms of "credible intervals". To this end 
we promote the model parameters to random variables and use the 
Bayes theorem to infer their posterior probability distribution of the 
parameters given the data and the model M : 



P(a \D,M) = 



C(D I a,M)P(a \ M) 
P(D j M) ' 



(19) 
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where £(D \ a, M) is the likelihood function, P(a j M) is the 
prior parameter model probability and P(D \ M) is the Bayesian 
"evidence". Notice that the likelihood is not a probability distribu- 
tion in the parameters a, since it gives the probability of the data 
D for a given value of the model parameters. In contrast, the ev- 
idence is the probability of the observed data within the assumed 
model M and appears as an overall normalization constant of the 
posterior. Hence, for parameter estimation purposes the evidence 
can be neglected since we are interested in finding the maximum 
of the posterior and the dispersion around it. This can be done by 
computing the likelihood function, which for a set of independent 
data reads as: 



[A ~ Dj(a) 
2 a? 



Central object of the Bayesian model selection is the evidence, 



C(D | a,M) = e~— = exp<^ - ^ 



(20) 

where Oi are the observational uncertainties and D t ! h (a) are the 
model predictions. However, notice that in our case, the position of 
the images, 9 a and 9b, play the role of independent variables in 
At ab and tab respectively. Since, the measured angles are known 
up to observational errors, it is important to propagate the angular 
uncertainties as well. This can be consistently done in the Bayesian 
framework. For a detailed discuss ion on the subjec t we refer to the 
exhaustive article by D'Agostini (D'A gostim1l2005h . 

Let assume that the observed image positions, 9 a and 9b, are 
determined with Gaussian uncertainties a a and as respectively. 
Let us indicate with &a and &b the true location of the images. 
Then, we can propagate the uncertainty on the independent vari- 
ables by marginalizing the likelihood over the angular error distri- 
bution. Assuming Gaussian errors this gives (for conciseness we 
drop the dependence on M): 

C(D | a) = jj f(& A ,& B )C(D | ® A ,&B,a)d® A d&B, 

(21) 



with 



/(©A,e B ) =exp 



E 



(gi - @, 

2a? 



(22) 



where in the integrand we have explicited the dependence of the 
likelihood C on the position of the images. If the angles are pre- 
cisely measured, then Eq. ((22} tends to a 5-Dirac and we recover 
the standard likelihood expression Eq. ( 120) . 



3.2 Bayesian model selection 

Parameter estimation only provides us with information on the 
quality of fit of a given model against the data. How do we choose 
between competing models characterized by different model pa- 
rameters? Preference solely based on the goodness-of-fit as mea- 
sured by the \ 2 f° r the best fitting model parameter values leaves 
us with partial information which is not consistently quantifiable 
in a probabilistic manner. Even then, how do we discount for the 
level of predictiveness of different models? More precisely, how 
do we decide whether data justify the choice of a more complex 
model Mi over one with a smaller number of free parameters Mo, 
while accounting for the extende d prior parameter space of M i ? 
As stressed in lLiddleetal.lj2007l) . such a question is related to the 
predictiveness of models rather than 'simplicity /complexity', since 
the former is not necessarily related to the number of free parame- 
ters. Parameter estimation does not address these questions, which 
require a further step already built in the Bayesian approach. 



P(D | M) 



da C{D | a, M)P(a | M), (23) 



which gives the probability of the data given the model M . Using 
the Bayes theorem we then have the probability of the model M 
given the data, 



P(M | D) oc P(D | M)P(M), 



(24) 



where P(M) is the prior belief on the model. Thus, given two 
competing models, Mo and Mi, we can base our preference on 
the ratio of the model probabilities given the same set of data 



P(Mo | D) _ P(D | Mo) P(Mo) 



P(Mi 
where the ratio 



D) 



Boi — 



P(D | Mi) P(Mi) 



P(D | Mo) 



P(D | Mi) 



(25) 



(26) 



is the so called "Bayes factor" of Mo to Mi. Supposing we have 
the same prior belief on the two models, then the Bayes factor gives 
us an estimate of the ratio of the model probabilities given the data. 
Since the evidence sets up a tension between the quality of fit of 
a given model and its prior predictiveness, then the Bayes factor 
accounts for the different size of the prior parameter space of the 
competing models. 

Bayes factors provide us with a self-consistent probabilistic 
measure to select a sample of lenses for which the observational 
data provide strong evidence in favor of a given model description. 
For instance, a simple lens model is the isothermal sphere which 
has no free parameters, Mo — {n = 2}. Using the Bayes factor we 
can compare it to a more complex description such as the power- 
law model with one free parameter Mi = {n}. In such a case 
we say that Mo is nested in Mi. Similarly, we can compare the 
power-law model with one free parameter Mo ~ {n} to the case in 
which a shear field is included with two additional free parameters 
Mi = {n,7, (f>j } . The additional presence of internal shear can be 
included in a third nested model M2 and confronted against Mo 
and Mi. For simplicity, we limit our current analysis to the latter 
models only. 

In order to assess the strength of evidence in favor of a model 
Mo, against a more complex one, Mi, w e use the Jeffreys scale 
Jjeffrevsll 196 lh as reported in Trotta (2007). To be as conservative 
as possible, we select lenses for which lnBoi > 5 and exclude 
those which provide strong evidence in favor of external shear, i.e. 
InBio > 5. We also exclude lenses for which we consider the 
model comparison to be inconclusive —5 < In Boi < 5. Notice 
that lnBoi = 5 corresponds to odds of 1 in about 150 in favor of 
Mo. The arrival of new data update these odds, further improving 
our knowledge of the lens system. 



4 TESTING BAYESIAN MODEL SELECTION 

In this section we perform a detailed analysis on a synthetic lens 
catalog to quantify the level of bias on the value of Ho inferred 
from a Bayes factor selected subsample of simulated lenses. To be 
as conservative as possible we consider lenses in the presence of 
external shear. This allows us to test to which extent Bayes factors 
can select those systems for which a simple power-law potential is 
sufficient to describe the simulated data and consequently evaluate 
the systematic bias on the inferred value of Hq. 
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Figure 1. Histogram showing the distribution of simulated values of 7 (left panel) and n (right panel) for the selected subsample So (solid red line) and the 
remaining dataset Si (green dash line). 




Figure 2. Histogram showing the distribution of values of Rab (left panel) and e (right panel) for the selected subsample So (solid red line) and the rejected 
dataset Si (green dash line). 



4.1 Synthetic Lens Catalogue 

We assume a flat ACDM model with fl m = 0.3 and a reduced 
Hubble consta nt h = 0.72. Using the publicly available software 
GRAVLENS jKeetodl201lh we generate a sample of 500 gravita- 
tional lenses at z; = 0.5 with sources at z s = 1.8. A more realistic 
simulation would be to distribute lenses and sources over a range 
of redshifts compatible with observations. However, this choice is 
not relevant to the purpose of this analysis and for simplicity we 
assume them to be at the same redshifts. 

The synthetic catalog is built by Monte Carlo generating a dis- 
tribution of values of n and 7 drawn from Gaussian distributions 
withn ~ Af(2, 0.3) and 7 ~ A/"(0.2, 0.1) consistently with current 
observational constraints discussed in Section 12.21 The position of 
sources is drawn from a uniform distribution in a cartesian angular 
map with coordinates — 1 < x, y < 1. For simplicity we fix b = 1. 
For each simulated lens, GRAVLENS computes the number, posi- 
tion and time-delay of the images. Out of 500 simulated lenses only 
312 have more than one image. In particular, 280 lenses with two or 
three images (double lenses) and 32 with four or five images (quad 
lenses). For simplicity we limit our analysis to the double lens sub- 
sample, since quad lenses require an accurate modeling of the lens 
angular structure which is essential to explain the formation of four 
images. We assume a 10% error on the time-delay, rather higher 



than the current mean error, and neglect angular errors on the im- 
age positions. 



4.2 Model Selection and Systematic Bias 

Using the data in the reduced catalog, we compute the likeli- 
hood function C over the prior lens parameter space of a simple 
power-law model A4o = {n} and one including external shear 
Mi = {n,7, </> 7 } respectively. We assume uniform priors with 
n G (1, 3) and 7 £ (0, 0.2). Then, we perform a numerical inte- 
gration to evaluate the bayesian evidence of each model and com- 
pute the Bayes factor, assuming a fixed cosmology with h = 0.7. 
Given the limited number of parameters we do not perform a Monte 
Carlo sampling of the likelihood, rather we compute it on a fine 
multi-dimensional grid. 

We select all lenses with mBni > 2.5. This gives us a sub- 
sample, So, consisting of 111 lenses for which image astrometry 
and time-delay provide strong evidence in favor of a simple power- 
law model description. We define the remaining lenses as the sub- 
sample Si. 

In Fig. Q] is shown the distribution of value of 7 (left panel) 
and n (right panel) for So (solid red line) and the remaining dataset 
Si (green dash line). As we can see the two subsamples cover the 
same range of values, thus indicating that Bayes factors are not sys- 
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Figure 3. Marginalized likelihoods on Ho from the analysis of So (solid 
red line), Si (green dot line) and the entire sample So + Si (dark blue 
dash line), all with model A^o- The dash-dotted light-blue line is obtained 
by analysing the entire sample with model M i . This model is undercon- 
strained by time-delay and astrometry data only, hence it is not surprising 
that due to the unbounded parameter degeneracies the inferred value of Ho 
strongly differs from the fiducial one. 



tematically selecting lenses with particular slope profiles or small 
external shear amplitude. This might seem at odd with using the se- 
lected subsample So to perform an unbiased parameter estimation 
using the simpler model A4q. However, as extens ively discussed in 
related literature (e.g. see Mukheriee et al.ll20060 . what the Bayes 
factors do is to put a tension between the capacity of the extended 
model A4 1 to better fit the data and its larger prior parameter vol- 
ume relative to the simpler model. Hence, even though the selected 
sample of lenses has non- vanishing shear, the data (image astrome- 
try and time-delay) do not justify the more complex model because 
the gain in fitting with external shear parameters is minimal com- 
pared to the size of the enlarged prior parameter space. The Bayes 
factors simply tell us that the simpler model A^o has to be pref- 
ered as it does a better job at describing the data with a smaller 
prior parameter space. By construction all simulated lenses have 
external shear and not surprising more than half of them are indeed 
discarded by the Bayes factor either as inconclusive or as favoring 
model M i . 

In Fig. [2] we plot the histogram of the values of Rab and e 
for So and Si respectively. We may notice that in So the distri- 
bution of values of Rab has a smaller scatter than for Si with 
Rab > 0.35. Similarly the distribution of values of e for the sub- 
sample So is narrower than Si and limited to values with e < 0.2. 
This clearly indicate that the Bayes factor analysis has selected 
lenses for which the angular configuration of the images is charac- 
teristics of double lens system whose astrometry an d time-delay s 
are less sensitive to external shear effect (e.g. see lQgurill2007h . 
hence successfully described by a simple model. Thus, the Bayes 
factors automatically performs an empirical model selection based 
on the structural properties of the lens systems. 

We are now in a position to quantify the bias on the inferred 
value of Ho from the selected lens subsample. Assuming model 
Mo we run a likelihood analysis on So, Si and their combination 
So + Si. The marginalized likelihoods on Ho are shown in Fig. [3] 
^From So we find Ho = 75.94jli 4 kms -1 Mpc -1 , quite remark- 
ably this recovers the fiducial value Ho = 72kms _1 Mpc _1 to 
within 3a statistical uncertainty of the best fit value, thus indicat- 
ing a systematic bias of ~ 5%. In contrast, the analysis of Si gives 
a highly biased value with Ho = 90. 45iJi km s~ 1 Mpc~ 1 , while 



the combined analysis gives Ho = 78.39 ± 1.5 kms _1 Mpc _1 . In 
principle, the residual ~ 5% systematic bias in So can be reduced 
or possible removed using additional lens measurements. Here, we 
have limited to use time-delay and astrometry data only to perform 
a proof of concept assuming simple lens toy models. However, in a 
real setup, the availability of the stellar kinematic measurements of 
the lens galaxy as well as high resolution imaging of the Einstein 
ring may provide additional constraints on the lens mass distribu- 
tion, thus reducing the internal lens model parameter degeneracy, 
which can lead to a more accurate lens model selection and conse- 
quently less biased results. In fact, the former provides an indepen- 
dent estimate of the lens mass at a radius different from that of the 
Einstein ring dTreu et al.f2 006). while the thickness of the Einstein 
ring i s sensitive probe of the density profiles (e.g. see ISuvu et all 
2009). This is to say that the residual systematic bias we have found 
from the analysis of the simulated lenses is not inevitable in strong 
lens cosmography. 

One may wonder whether analysing the entire dataset using 
model A^i may lead to less biased results given that the data 
have been generated using such a model. However, we should re- 
mind that model Mi is underconstrained using double lenses. This 
can introduce very large parameter degeneracies and eventually 
strongly bias the cosmological parameter inference. To this puipose 
we have run a likelihood analysis of the whole sample So + Si 
assuming model Mi. The marginalized likelihood on Ho is plot- 
ted in Fig. [3] (dash-dotted light-blue line) from which we obtain 
Ho = 44 ± 1.5 kms _1 Mpc _1 . This shows how subtle the param- 
eter inference can be. In a real dataset we will not know a priori 
whether external shear is indeed present or not. Thus, given a set 
of data which we aim to use to extract cosmological parameter in- 
formation we are better guided by Bayesian model selection than 
blindly adding parameters to model the complexity behind the data. 
However, as already stressed above we should not dispair since ad- 
ditional constraints from independent mass proxies may help re- 
duce or break internal model parameter degeneracies and conse- 
quently the robustness of the lens model selection. 



5 APPLICATION TO REAL LENS DATA 
5.1 Data sample 

Gravitational time-delays have been measured in 21 strong lens 
systems. Out of this sample we only consider double image lenses 
in which a far distant quasar is lensed into two images by a fore- 
ground galaxy. This reduces our initial dataset to 12 lenses, whose 
characteristics are quoted in TableQ] The object responsible for the 
lensing has been unambigously detected in all listed lenses. How- 
ever, current observations do not provide clear indications whether 
the lensing galaxies are part of a group/cluster or whether pertur- 
bators are present along the line of sight. Thus, the presence of 
external shear cannot be a priori excluded, potentially leading to 
bias selection effects. 

For each lens in the sample we compute the asymmetry Rab , 
the angular separation Oab and the deviation from colinearity e 
as indicators of perturbations about a monopole-like lens potential. 
These are also given in Table Q] while in Table IbTI in Appendix |B1 
we report the astrometry. 
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Table 1. Double-image lenses. z; and z s are the lens and source redshift respectively, 8i are the angular position of the images relative to the lens position in 
arcsec, At is the time-delay (in days) and F^/Fj the flux-ratios. Observational uncertainties are 1 a errors. We also quote the derived values of the asymmetry 
Rab, angular separation ©ab and deviation from colinearity e. 



Lens 


pair (ij) 


z\ 




2s 


6i (") 


8j (") 


At = ti - ij 


F i 


Rab 


©AS 


e 


B 1600+434^ 


dA 


U.OOO 


U 




U.Uo I ± U.UU4 


U.zoU ± U.UUo 


+ 1U.0 ± U.z 


U.zo ± U.UUo 


U.oo 


ZUj.O 


n 1 A 
U. i4 


AB 


0.414 


1 


589 


1.14 ±0.075 


0.25 ±0.074 


-51.0 ±2.0 


1.75 ±0.34 


0.45 


200.2 


0.11 


FBQ 0951+2635^ 


AB 


0.260 


1 


246 


0.886 ± 0.004 


0.228 ± 0.008 


-16.0 ±2.0 


3.15 ±0.05 


0.59 


201.2 


0.12 


HE 1104-1805£ 


AB 


0.729 


2 


319 


1.099 ± 0.004 


2.095 ± 0.008 


±152.2 ± 3.0 


2.84 ±0.06 


0.32 


175.9 


0.02 


HE 2149-2745^ 
PKS 1830-21 1£ 


AB 


0.603 


2 


033 


1.354 ±0.008 


0.344 ±0.012 


-103.0 ± 12.0 


4.0 ±0.5 


0.59 


178.9 


0.006 


AB 


0.89 


2 


507 


0.67 ±0.08 


0.32 ±0.08 


-26 ±5 


1.03 ±0.02 


0.36 


199.5 


0.11 


Q0142-100S. " 


AB 


0.49 


2 


719 


1.855 ±0.002 


0.383 ± 0.005 


-89 ±11 


6.3 ±0.1 


0.66 


167.8 


0.07 


00957+561^ 


AB 


0.36 


1 


413 


5.220 ±0.006 


1.036 ±0.11 


-417.09 ± 0.07 


1.35 ±0.04 


0.67 


154.5 


0.14 


SBS 0909+532^, 


AB 


0.830 


1 


377 


0.415 ±0.126 


0.756 ±0.152 


±45.0 ±5.5 


0.32 ± 0.03 


0.29 


220.3 


0.22 


SBS 1520+53u£ 


AB 


0.717 


1 


855 


1.207 ±0.004 


0.386 ± 0.008 


-130.0 ±3.0 


2.9 ±0.4 


0.52 


202.6 


0.13 


SDSS 11206+4332^ 
SDSS J 1650+425 1£ 


AB 


0.748 


1 


789 


1.870 ±0.088 


1.278 ±0.097 


-116 ±5 


0.74 ±0.2 


0.19 


132.9 


0.26 


AB 


0.577 


1 


547 


0.872 ±0.027 


0.357 ±0.042 


-49.5 ± 1.9 


6.2 ±0.31 


0.42 


145.8 


0.19 



a Image pos i tions fro m iLehar et aT] (2000) computed relative to the l ens center as mea sured by lYork et all (|2005|); lens and source redshift s from 
iBrowne et alll 19931) and Cohen et aU 12003l) respectively, time-delay from Biggs e t al.1 J 1999b and flux ratio from lWucknitz. B iggs & Browne 1 20041) . 
" Image positions and redshifts from Koop mans. de Bruvn & Jacksonl jl998l) . time- delay fromlBurud et al I j2000t) . flux-rati o fromlDai & Kochanek l2005t) . 
c Image positions and redshifts as listed in iKochanek et alj <2008l) . time-delay from Jakobsso n et al.l J2005h - flux-ratio from Shalyapin et al. (2009). 
^ Image positions and redshifts from Lehar et al. ( 2000 ), time-del ay and flux-ratio from lPoindexter et al.ld2007| 



Image positions and redshifts as listed in Kochanek et al. 1 2008), time-delay and flux-ratio from Burud 



(2007), 

et al.N2002al) . 



f Image positions and redshifts from Meylan et al. (2005), time-delay from lLovell et al J j 1 998h - fl ux-ratio fro m Courbin et alj fl9 
9 Image positions and redshifts as listed in lKochanek et alj 120081) . time-delay and flux-ratio from Koptelova et all l2012l) l 

Image positions and redshifts as listed in Kochanek et al. (2008), time-delay from lCollev et alj d2003l) . flux-ratio from lHaarsma et al II I 19991) . 
* Image positions and redshifts as listed in Kochanek et al. 12008), time-delay from Ullan et al. (2006), flux-ratio from Dai & Kochanek 1 2009). 
3 Image positions and redshifts as listed in Koch anek et alj (2008), time-delay from lBurud et alj (2002b ). flux-ratio from Auger et al. (2008). 
™ Image positions and redshifts as listed in Oguri et al (2005), time-delay and flux-ratio from Paraficz, Hiorth & Eh'asdottir (2009). 
' Image positions and redshifts as listed in Kochanek et al. ( 2008), time-delay and flux-ratio from Vuissoz et al. ( 2007). 



5.2 Model parameter priors 

We want to test whether the data listed in Table[T]provide strong ev- 
idence for or against a simple power-law model, Mo = { n }, or one 
which includes external shear, Mi = {n, 7, 4>~t}- Uninformative 
data cause the Bayes factors to depend on the size of the assumed 
prior parameter space. To test for such dependence we distinguish 
two different sets of priors for each lens model parameter. 

To be as conservative as possible we assume a uniform "large 
prior" on the radial slope of the potential, n £ (0, 3) and a uni- 
form "small prior" corresponding to 11 £ (1,3). The former in- 
cludes shallow galax y profiles (n < 1) such as those found in e.g. 
Salucci et aD {2007), while the latter is limited to the cuspy pro- 
files, including the isothermal sphere, which are typical of Cold 
Dark Matter halos. Similarly, we assume flat priors on 7 and con- 
sider two prior intervals, 7 £ (0, 0.1) and 7 £ (0, 0.2), while the 
external shear orientation angle is assumed </> 7 £ (0, n). 

In the evaluation of the Bayes factors we consider a flat con- 
cordance ACDM model and assume for si mplicity a value com- 
patible with WMAP-7yr dLarson etal.l201lh on the matter density, 
f^m = 0.27. In order to test the dependence of the Bayes factors 
on the value of the Hubble constant we assume three prior values: 
H = 65, 70.2 and 75 kms^Mpc" 1 respectively. 

5.3 Bayes factor analysis 

The results of the numerical computation of the Bayes factors ob- 
tained from the analysis of time-delays and image astrometry are 



summarized in Fig. [4] while in Fig. [5] we show the results obtained 
including flux-ratio measurements. For each lens in the sample we 
mark the value of log 10 B01 under different lens model parameter 
priors and for values of the Hubble constant Ho = 65 (left panel), 
70.2 (middle panel) and 75kms _1 Mpc _1 (right panel) respec- 
tively. The hatched areas correspond to values of the Bayes factors 
for which the data provide strong evidence in favor of the simple 
power-law description (In Box > 5), the presence of external shear 
(In Bio > 5). We assume the Bayes model comparison to be incon- 
clusive in the region — 5 < In B01 < 5 which is a very conservative 
cut. 

Let us focus on Fig. [4] First we can see that Bayes factors 
show no variation for the large and small priors on n, except 
for Q0957+561. On the other hand, most lenses exhibit prior de- 
pendence on 7, which can vary in strenght with the cosmology. 
FBQ 0951+2635, HE 2149-2745 and Q0957+561 are the only ones 
for which the prior on 7 seems to have no influence whatsoever. 
We may also notice a weak dependence of the Bayes factors on 
the value of h. In particular, B1600+434 and SDSS J1650+4251 
have Bayes factor which are clearly above our conservative cut for 
h = 0.702 and h = 0.75, while these shift in the upper bound on 
the inconclusive interval for h = 0.65. However, since the value 
of Bqi is still rather high, these lenses would pass a slightly less 
conservative cut even in the worst case scenario of h = 0.65. 

All lenses except HE 1104-1805 and Q0957+561 have Bayes 
factors favoring Mo, though they may lie in the inconclusive area. 
Q0957+561 is the only lens clearly favoring model Mi in some 
case. One final remark concerns SBS 1520+530. Our analysis in- 
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Figure 4. Bayes factors for the sample of lenses listed in Table[T]and obtained from the analysis of the image positions and time-delays for different lens 
model parameter priors and assuming a hard prior on the reduced Hubble parameter with h = 0.65 (left panel), h = 0.702 (middle panel) and h = 0.75 
(right panel). 



Iff 



10' 



10" 



10 : 



-10'" 



-10' 



-Iff 







h = 0.75K N 








U 

















Prior on y: <0.1 <0.2 <0.1 <0.2 
Prior on n: large small 



<0.1 <0.2 <0.1 <0.2 
large small 



<0.1 <0.2 <0.1 <0.2 
large small 



B02 18+357 
B 1600+434 
FBQ 095 1+2635 
HE 1104-1805 
HE 2149-2745 
PKS 1830-211 
Q0 142- 1 00 
00957+561 
SBS 0909+532 
SBS 1520+530 
SDSS J 1206+4332 
SDSS J 1650+4251 



Simple model 
favored 




Model with 
shear favored 



Figure 5. As in Figure |4]including information from flux-ratio measurements. 



dicates that the image position and the time-delay of this lens 
strongly favor a simple power-law lens model without external 
shear. This is not in contr ast with observations of the lens proper- 
ties bv lAuger et al.l d2008h who have shown the presence of galaxy 
groups at nearby redshifts on the lens line of sight. It is plausi- 
ble that mass potential is centered on the lens (as suggested by the 
galaxy-groups position relative to the lens) thus implying a negli- 
gible quadrupole contribution. Hence, as discussed in section |4]the 
extra parameters from modeling the external shear may provide a 
minimal gain in fitting the data and in such a case the Bayes factor 
favour s the simpler model description. Furthermore, Au ger et all 
( 2008) find that the system can be well fitted by a nearly isothermal 
density profile. 

In Fig.[5]we plot the Bayes factors obtained from time-delays 
and flux-ratio measurements. As we can see including flux-ratios 
shifts the Bayes factors of several lens systems in the range which 



strongly favor the presence of external shear. However, it is possible 
that the extra parameters of model M\ fit unmodeled systematics 
affecting the flux-ratios rather than the effects of an external shear 
field. As already mentioned in Section l2~2l current observations al- 
ready indicate that dust extinction as well as microlensing and sub- 
structure alter the flux-ratios. Since we are far from being able to 
quantitatively account for such contaminations, it is premature to 
infer any conclusion based on these data. 

In the light of the analysis of Bayes factors inferred from po- 
sition of the lens images and time-delays we select B 1600+434, 
SBS 1520+530 and SDSS J1650+4251 out of the initial lens sam- 
ple. These are the only lenses which pass our selection criterion 
independently of the model parameter priors with the power- law 
model strongly favored over that including external shear with odds 
of more than 150 to 1. In the next Section we will present the con- 
straints on Ho obtained from the likelihood analysis of these lenses. 
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Figure 6. Marginalized 1 and 2<r contours in the Q m — 71 plane without 
WMAP prior for B 1600-434. 



Fi gure 8. din | H m , i^o)/£max for each of the selected lenses. 




Figure 7. As in Fig.|6]for the Ho — n plane without WMAP prior. 



6 COSMOLOGICAL PARAMETER INFERENCE 

Using the image positions and time-delays of B 1600+434, 
SBS 1520+530 and SDSS 11650+4251, we perform a likelihood 
analysis to infer constraints on Ho after marginalizing over the 
angular errors of the image locations, the slope of the power-law 
model of each lens m and the matter density fi TO . We assume flat 
priors on 1 < n < 3 and 30 < H a [kms _1 Mpc _1 ] < 150. 
For each lens we compute the likelihood function C(rn, fl m , Ho) 
as given by Eq. d21t . We illustrate the parameter degeneracies for 
B 1600+434, in Fig.|6]and Fig.[7]we plot the 1 and 2a marginalized 
contours in the planes fi m — n and Ho — n respectively. We can 
see that there is no correlation between n and fl m for such dataset, 
while a more pronounced degeneracy i s present between n and Ho, 
as was already pointed in lSuvul (|2pi2). Introducing a prior on fl m 
does not change the likelihood contour significantly, as there is lit- 
tle degeneracy between fi m and either n or Ho- 

In Fig. [8] we plot the 1 -dimensional marginalized likelihoods 
in the prior interval of n for each of the lenses obtained assuming a 
Gaussian p rior on O m = 0.26 6 ± 0.029 consistent with WMAP7- 
yrs results dLarson et al.ll201 if) . We can see that data can only pro- 
vide a lower limit on the value of the slope of the lens potential, 
n > 1.5. Similarly, in Fig. |9]we plot the marginalized likelihoods 
of Ho for each lens. We may notice that the likelihoods peak around 
nearly the same value of Ho and have a large dispersion around the 
maximum value. 

The combined constraints on Ho are inferred by computing 
the total likelihood sample: 

C(Q m , Ho) = Yl [ C{m, Am, H ) P(m) dn,, (27) 

where the sum is over the likelihood of each lens weighted by 
P(rii) the uniform prior in the interval 1 < m < 3. Hence, we 
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Fi gure 9. Ci^Ho I ^m, ^)/^max for each of the selected lenses. 



do not assume a hard prior on n at the best-fit value of each lens, 
rather we impose the model and propagate its parameter uncertain- 
ties by marginalising the likelihood over them. In Fig.[l0]we plot 
the 1 and 2a contours in the fi m — Ho plane (left panel) and the 1- 
dimensional marginalized likelihood for Ho (right panel) with and 
without WMAP7-yrs prior. As we can see, time-delays are mostly 
insensitive to Q m . The average and standard deviation values are 



Ho 



Tei^kms^Mpc- 1 and Ho 



78] 



''kms -1 Mpc _1 



without and with WMAP7-yr prior respectively. Notice that the 
marginalized likelihood is non-Gaussian and strongly skewed to- 
wards small values of Ho - Thus, the average values differ the max- 
imum ones which are at Ho = 85 km s~ Mpc -1 and Ho = 
88kms _1 Mpc -1 without and with WMAP7-yr prior respectively. 

It is worth highlighting that, while on the observational sam- 
ple, about a quarter (three out of twelve) of the lenses are selected 
through our analysis, in the simulated sample this ratio is higher 
than a third (111 out of 280), which seems to indicate that observa- 
tional samples are more contaminated than our simulated sample, 
leading to possibly higher bias if the whole sample is used to mea- 
sure Ho- 



7 CONCLUSIONS 

In this study we have used Bayesian model selection techniques to 
determine the lens mass model which given the data has the high- 
est probability to describe a strong gravitational lens system. Rather 
than modeling a lens in all its complexity, our goal is to focus on 
selecting a model whose parameters can significantly influence the 
time-delay such as to infer unbiased cosmological constraints when 
averaging individual mass model parameter uncertainties on a ho- 
mogeneous lens sample. To this end we have focused on double 
lenses and used Bayes factor statistics to select a sample of lenses 
that can be described by a power-law model without external shear. 
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which the quadrupole structure induces an internal shear that needs 
to be taken into account to correctly describe the lens time-delay . 
The lens reconst ruction technique described in [Alard (2007 . [200a) 
and adapted in Habara & Yamamoto ( 1201 ll) could be particularly 
useful in this case especially in combination with the Bayesian se- 
lection approach described here to determine the order of the per- 
turbative lens reconstruction. Though, such technique may present 
difficulties since in most cases we do not see arcs but only point- 
like images. 

The number of known lensed quasars which are potentially 
good targets for time-delay measurements is currently of order a 
hundred. In the future such dataset will increase thanks to numer- 
ous survey programs. The use of the Bayesian model selection 
we have discussed can provide large homogeneous subsample of 
lenses that are suitable for unbiased c osmological paramet e r infer - 
ence. In particular, along the line of IParkinson & Liddla d2010h . 
using Markov Monte Carlo Chain techniques one could perform 
a Bayesian model averaging of cosmological parameters over ho- 
mogeneous subsamples of lenses each described by a class of lens 
models. 



Figure 10. Top panel: 1 and 2cr contours Q m — Ho plane from the com- 
bined analysis of the selected lens sample. Bottom panel: C(Hq) / Cmux 
for the same sample. 

We have tested this approach on a simulated sample of 500 lenses. 
The likelihood analysis of the selected subsample assuming no ex- 
ternal shear recovers the fiducial value of the Hubble constant with 
3a statistical uncertainty, thus a systematic selection bias is no 
greater than ~ 5%. On the other hand, since the more complex 
model including external shear parameter is underconstrained us- 
ing double lenses, the likelihood analysis result in strongly biased 
results due to the marginalized effects of large lens model parame- 
ter degeneracies. Therefore, a given set of data should not be brute 
force analysed assuming the model capable of accounting for the 
maximal complexity. Rather, performing a Bayes factor analysis 
may provide a more effective guide to data model selection. 

The application to a sample of observed lenses indicates that 
out of the initial dataset, nine have Bayes factors favoring a sim- 
ple power-law lens model, though six of them lies in the "in- 
conclusive" interval and do not pass our conservative cut. These 
are B0218+357, FBQ 0951+2635, HE 2149-2745, PKS 1830-211, 
Q0142-100, and SBS 0909+532. It is possible that more accurate 
time-delay measurements will update their Bayes factor value and 
give us a better knowledge of the appropriate lens model descrip- 
tion. We have performed the same analysis including information 
from lens flux-ratio measurements. In such a case, the Bayes fac- 
tors indicate that external shear must be taken into account for 
most of the lenses. However, several source of astrophysical sys- 
tematic errors may affect the flux-ratios such as dust extinction and 
microlensing events that do not reflect the lensing magnification 
caused by the lens mass distribution responsible for the images and 
time-delays. Indeed, uncontaminated flux-ratios will provide inter- 
esting additional constraint on the lens model. Nevertheless, it is 
worth reminding that time-delays are less sensitive to substructure 
than flux-ratios and require a less complex modelling. 

Here, we have restricted the analysis to double lens systems. 
Such systems have fewer constraints than quads, nevertheless the 
advantage is that they can be described in terms of simpler mod- 
els with fewer parameters. This is not the case for quad lenses in 
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APPENDIX A: MODEL WITH EXTERNAL SHEAR 

In this Appendix we present a derivation of the time-delay and mag- 
nification equations respectively in the case of a power-law poten- 
tial with external shear. The latter reads as 



7y cos2( 



r)- 



(Al) 



We can relate one of the model parameters to the remaining ones 
using the lens equation, since both images result of the same source 
at (3. This gives us the following equation: 



I/3(0b)| 



(A2) 



which can be rewritten, using the lens equation (3 = 8 — \jif)(6), 
as a second order equation in X = & n_1 : 



with 

U 



VX 1 + 2VX + W = 



g 2(2-n) 



V = ^-' 1 (1 + 7Cb)-^"(1 + 7 Ca) 
W = 0a(1 + 27Ca + J 2 ) - #1(1 + 27Cb + 7 2 ) 



(A3) 

(A4) 
(A5) 
(A6) 
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where Ca.b = cos2(<^>a,b — <t>~t)- Assuming that 8a > 9b and 
n > 1, and using the usual notation A = V 2 — UW we have the 
following solutions: for A > 



Table Bl. Astrometry for all lenses used in the data analysis 



the solution b n 



-v + 7a 
u 



(A7) 



is eliminated since it does not reduce 



(A8) 



to the right value for 7 = (note however that it would be the 
correct solution if n < 1). Note also that for n = 1, there is no 
parameter b and this equation has no meaning; for A = 0, the 
solution is simply 

-W 

IV '■ 

while for A < 0, there are no real solutions and the lens cannot be 
described for this combination of n, 7 and 7 . Having reduced the 
number of lens model parameters the time delay between images A 
and B can be computed using the value of tx and ta I 



U = (1 + zi) 



~ 1 V ^ s 2 +m: 



(A9) 



where we have used the lens equation /3 — 6>i = — v ^(fii). In the 
case of the magnification we have for each image: 



1 - (2 - n) 



+(1 - n)yQ [ 



(A10) 
(All) 



and the flux ratio is simply given by Fab = /xa//xb- 



APPENDIX B: LENS DATA ASTROMETRY 

In table lBll we report the astrometric data of the positions of images 
(A,B) and lens (L) for all lenses used in the lens data analysis. The 
references to the data are listed in table Q] 



lens 


obj ect 


Aft 


A<5 


B02 18+357 


L 


EE 


EE 




A 


-0.250 ± 0.005 


-0.125 ±0.007 




B 


0.057 ±0.006 


0.001 ±0.008 


B 1600+434 


L 


ee 0±0.05 


ee 0±0.05 




A 


-0.33 ±0.01 


1.09 ±0.01 




B 


0.07 ±0.01 


-0.24 ±0.01 


FBQ 0951+2635 


G 


0.760 ± 0.003 


-0.455 ±0.003 




A 


EE 


EE 




B 


0.900 ± 0.003 


-0.635 ± 0.003 


HE 1104-1805 


G 


0.974 ± 0.003 


-0.510 ±0.004 




A 


EE 


EE 




B 


2.901 ±0.003 


-1.332 ±0.003 


HE 2149-2745 


G 


0.714 ±0.007 


1.150 ±0.005 




A 


EE 


EE 




B 


0.890 ± 0.003 


1.446 ± 0.003 


PKS 1830-211 


G 


0.498 ± 0.004 


-0.456 ± 0.004 




A 


EE 


EE 




B 


0.649 ± 0.001 


-0.724 ±0.001 


Q0142-100 


G 


1.764 ±0.003 


-0.574 ±0.003 




A 


EE 


EE 




B 


2.145 ±0.003 


-0.613 ±0.003 


Q0957+561 


G 


1.406 ±0.006 


-5.027 ±0.005 




A 


EE 


EE 




B 


1.229 ±0.005 


-6.048 ± 0.004 


SBS 0909+532 


G 


0.415 ±0.125 


-0.004 ± 0.081 




A 


EE 


EE 




B 


0.987 ±0.003 


-0.498 ± 0.003 


SBS 1520+530 


G 


1.141 ±0.003 


-0.395 ± 0.003 




A 


EE 


EE 




B 


1.429 ±0.003 


-0.652 ± 0.003 


SDSS J1206+4332 


G 


-0.664 ±0.137 


1.748 ±0.028 




A 


EE 0± 0.011 


EE 0± 0.010 




B 


-0.098 ± 0.006 


2.894 ± 0.009 


SDSS J1650+4251 


G 


0.017 ±0.032 


-0.872 ±0.026 




A 


EE 


EE 




B 


0.223 ±0.002 


-1.163 ±0.001 
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